1 Data set up

load(file.path(data_path, "Candies.RData"))

pander("table(design$Judges, design$Candies)")

table(design\(Judges, design\)Candies)

pander(table(design$Judges, design$Candies))
  1 2 3 4 5
01 3 3 3 3 3
02 3 3 3 3 3
03 3 3 3 3 3
04 3 3 3 3 3
05 3 3 3 3 3
06 3 3 3 3 3
07 3 3 3 3 3
08 3 3 3 3 3
09 3 3 3 3 3
10 3 3 3 3 3
11 3 3 3 3 3
raw_outcomes <- outcomes



df <- data.frame(raw_outcomes, design)
library(doBy)



# Mean attributes values by judges
table_summary_mean <- summaryBy(. ~ Judges , data = df, 
  FUN = mean )

pander(table_summary_mean)
Table continues below
Judges Transp.mean Acid.mean Sweet.mean Raspb.mean Sugar.mean
01 9.29 4.69 7.83 4.43 5.75
02 8.36 3.69 4.2 4.49 5.73
03 9.95 4.78 6.99 7.07 5.7
04 9 3.22 4.28 3.73 6.1
05 5.96 7.59 9.27 7.92 4.72
06 8.67 5.55 7.02 5.73 4.81
07 7.64 10.36 7.68 5.99 5.19
08 8.26 4.17 6.96 4.2 5.67
09 8.29 5.79 7.69 6.07 6.04
10 8.76 4.92 3.98 5.11 5.98
11 8.35 3.39 7.77 6 5.31
Bites.mean Hard.mean Elastic.mean Sticky.mean
10.08 9.21 8.54 8.87
6.48 7.14 7.93 7.33
9.12 8.83 7.53 7.55
9.66 9.57 9.32 7.26
8.5 8.96 9.02 7.08
9.13 9.46 8.79 8.89
9.39 9.26 9.45 9.49
8.88 8.36 7.75 8.79
9.7 9.35 9.1 9.41
9.5 8.68 9.68 8.29
9.58 9.53 8.58 7.79
res <-as.matrix(table_summary_mean[,-1])

# par(xpd=TRUE, mar =c(2,2,2,6))
# barplot(res, beside = TRUE, col=rainbow(11), 
#         main = "mean rating per judge and attribute")
# legend("topright", inset = c(-0.2,0), legend = paste0("As", 1:11), 
#        col = rainbow(11), pch=15)

2 remove outliers ??

exclude <- FALSE

if(exclude){
  index <- which(!rownames(raw_outcomes) %in% c("0211", "0713"))
  raw_outcomes <- raw_outcomes[index,]
  design <- design[index,]
}
Candies <- design$Candies
Judges <- design$Judges

2.1 Graph

# pdf(file.path(fig_path, "SDA_outcomes.pdf"), width = 8, height = 5,
#     pointsize=19)
par(mar=c(4,4,2,5), xpd = TRUE)
col <- gg_color_hue(n=5)[design$Candies]
outcomes <- as.matrix(outcomes)
plot( outcomes[1,], type="l",  xaxt="n",
     ylim=range(outcomes), col=col[1],
     xlab="Attribute", ylab = "Rating", main = "Sensory Data outcomes")
for (i in 2:dim(outcomes)[1]){
  lines(outcomes[i,], col=col[i])
}
axis(side=1, at = 1:9, labels = colnames(outcomes))
legend("topright", legend=levels(design$Candies), 
       col=rainbow(n=5), lwd=1, title="Candies", 
       inset=c(-0.24,0), box.col = "white")

# dev.off()

2.2 Hotelling T^2 for outlier detection

model = pca(raw_outcomes, scale = FALSE, 
            info = 'Simple PCA model', lim.type = "jm")
# model$ncomp.selected
ncomp <- 4

# plotVariance(model)

# plotResiduals(model, show.labels = TRUE, ncomp = ncomp, 
#               main = "Squared residual distance vs Hotelling T2 distance", 
#               norm = TRUE)

# plotResiduals(model, show.labels = TRUE, ncomp = ncomp, 
#               main = "Squared residual distance vs Hotelling T2 distance")

Qlim <- model$Qlim
T2lim <- model$T2lim
rownames(Qlim)[2] <- rownames(T2lim)[2] <- "Out_limit"



# In case of PCA the critical limits are just shown 
# on residual plot as lines and can be used for detection 
# of extreme objects (solid line) and outliers (dashed line).

plot_hotelling <- function(){
  xlim <- range(model$calres$T2[,ncomp])
  xlim[1] <- xlim[1]*0.9
  xlim[2] <- xlim[2]*1.1
  ylim <- range(model$calres$Q[,ncomp])
  ylim[1] <- ylim[1]*0.9
  ylim[2] <- ylim[2]*1.1
  plot(model$calres$T2[,ncomp], model$calres$Q[,ncomp], 
     main = "Diagnostic plot for score \nand residual outliers", xlab = "Hotelling T2 distance", 
     ylab ="Squared residual distance", pch = 16, xlim = xlim, 
     ylim = ylim, cex.lab=1.25, cex.main=1.3)
  abline(h=Qlim["Out_limit",ncomp], v=T2lim["Out_limit",ncomp], lty =3)
  legend("topright", legend = "Outlier limit", lty = 3, cex=1.2)
  index1 <- which(model$calres$T2[,ncomp]>=T2lim["Out_limit",ncomp])
  index2 <- which(model$calres$Q[,ncomp]>=Qlim["Out_limit",ncomp])
  index_ho <- unique(c(index1,index2))
  text(x = model$calres$T2[index_ho,ncomp], y=model$calres$Q[index_ho,ncomp], 
       labels = names(model$calres$T2[index_ho,ncomp]), pos = c(1,2,3,4))
}

plot_hotelling()

# pdf(file.path(fig_path,"SDA_hotelling.pdf"), height = 6, width = 6,
#     pointsize = 12)
plot_hotelling()

# dev.off()

3 Distribution of the raw outcome variables

par(mfrow=c(1,3))
for (i in 1:ncol(raw_outcomes)){
  # qqPlot(raw_outcomes[,i])
  d = density(raw_outcomes[,i])
  # hist(raw_outcomes[,i])
  plot(d)
}

4 Dimension reduction by PCA

n = nrow(raw_outcomes)

###################################################
# Dimension reduction by PCA
###################################################

# ===== PCA ===== #

res_pca <- MBXUCL::SVDforPCA(x = raw_outcomes)

pander("Cumulated variance")

Cumulated variance

pander(res_pca$cumvar)
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
74.69 84.08 90.55 93.44 95.42 97.3 98.58 99.67 100
barplot(res_pca$var, main="screeplot", ylab="% var", ylim=c(0,80))

4.1 PCA plots

df <- data.frame(PC = as.character(1:length(res_pca$var)),
                 var = res_pca$var)

screePlot <- ggplot(df, aes(y=0,yend=var,x=PC,
                          xend=PC))+ geom_segment()+
          labs(title= "Scree \nplot",
                        x = "PC", y="% var")+theme_classic()+
  theme(axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12))

index <- 27
pch <- c("1"=15, "2"=2, "3"=19, "4"=4, "5"=8)

scoresPlot <- DrawScores(res_pca, axes=c(1,2),drawNames = FALSE, 
           main  = "Scores plot", 
           pch = Candies, size=2, 
           legend_shape_manual  = pch) +
   coord_cartesian(xlim = c(-25, 25), ylim = c(-15, 15)) +
  annotate("text", y = (res_pca$scores[index,2] +
                                1*c(1)),
                  x = res_pca$scores[index,1],
                  label = rownames(res_pca$scores[,1:2])[index])
## Scale for 'shape' is already present. Adding another scale for 'shape', which
## will replace the existing scale.
loadPlot <- ScatterPlot(res_pca$loadings[,1], res_pca$loadings[,2], 
                        points_labs = rownames(res_pca$loadings), cex.lab=4) +
  geom_vline(xintercept = 0, lwd=0.1) + 
  geom_hline(yintercept = 0, lwd=0.1) + 
  xlab(label=paste0("PC1 (", round(res_pca$var[1],2),"%)")) +
         ylab(paste0("PC2 (", round(res_pca$var[2],2),"%)")) +
   coord_cartesian(xlim = c(-0.55, 0.55), ylim = c(-0.8, 0.8)) +
  ggtitle(label = "Loadings plot")

4.2 Apply first step of LiMMPCA (PCA)

res_pca <- dimreducPCA(data = raw_outcomes, pcvar = 99)

nPC <-dim(res_pca$pca_scores)[2]
spectra_PCA_scores <- res_pca$pca_scores
spectra_PCA_loadings <- res_pca$pca_loadings

centeringFactor <- colMeans(outcomes)

# outcomes
outcomes <- spectra_PCA_scores
rownames(outcomes) <- rownames(design)

4.3 Graph

p <- ggarrange(screePlot, loadPlot,
                ncol = 2, nrow = 1,
                widths  = c(0.3,0.95))

PCA_plots <- ggarrange(p, scoresPlot,
                ncol = 1, nrow = 2)

PCA_plots

# ggsave(file.path(fig_path,"SDA_PCA_outcomes.pdf"), plot = PCA_plots,
#          height = 10, width = 6.5, scale=0.75)

5 Tests with univariate ANOVA

mat <- cbind(subject=1:dim(design)[1], design,  outcomes)
mat$subject <- as.factor(mat$subject)

########################
# for one response 
########################

# res_aov <- aov(PC1 ~ Candies*Judges + Error(subject/Judges), data=mat)

# sumar <- summary(res_aov)

# sumar <- sumar[[1]][[1]]

#                   Df  Sum Sq Mean Sq  F value  Pr(>F)    
# Candies             4 31470.6  7867.7 780.1762 < 2e-16 ***
# Judges          10   224.9    22.5   2.2304 0.02089 *  
# Candies:Judges  40   707.7    17.7   1.7545 0.01158 *  
# Residuals         110  1109.3    10.1                     
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


# MeanSq <- sumar$`Mean Sq`
# names(MeanSq) <- gsub(" ", "", rownames(sumar))

# MeanSq["Candies"]
# MeanSq["Judges"]
# MeanSq["Candies:Judges"]
# MeanSq["Residuals"]
# 
# MeanSq["Candies"]/MeanSq["Candies:Judges"]
# MeanSq["Judges"]/MeanSq["Residuals"]
# MeanSq["Candies:Judges"]/MeanSq["Residuals"]

########################
# for all the responses
########################
sumar_list <- vector(mode="list", length=nPC)

pvalues <- matrix(data = NA, nrow=nPC, ncol=3, 
                  dimnames = list(NULL, c( "Candies","Judges","Candies:Judges")))

for (i in 1:nPC){
  form <- paste(paste0("PC",i), "~ Candies*Judges + Error(subject/Judges)")
  res_aov <- suppressWarnings(aov(as.formula(form), data=mat))
  sumar_list[[i]] <- summary(res_aov)[[1]][[1]]
  cat("PC ", i)
  print(sumar_list[[i]])
 
 # pvalues
  nam <- rownames(sumar_list[[i]])
  nam <-  trimws(nam , which = c("both", "left", "right"))
  meanSq <- sumar_list[[i]]$`Mean Sq`
  names(meanSq) <- nam
  Df <- sumar_list[[i]]$Df
  names(Df) <- nam
  
  pval_Candies <- pf(meanSq["Candies"]/meanSq["Candies:Judges"],
                     df1=Df["Candies"], df2=Df["Candies:Judges"], 
                     lower.tail = FALSE)
  pval_Judges <-  pf(meanSq["Judges"]/meanSq["Residuals"], 
                     df1=Df["Judges"], df2=Df["Residuals"], 
                     lower.tail = FALSE)
  pval_CA  <-  pf(meanSq["Candies:Judges"]/meanSq["Residuals"],
                  df1=Df["Candies:Judges"], df2=Df["Residuals"], 
                     lower.tail = FALSE)
  pval <- c(Candies=pval_Candies,
               Judges=pval_Judges, 
               CA=pval_CA)
  pvalues[i,] <- pval

}
## PC  1                Df  Sum Sq Mean Sq  F value  Pr(>F)    
## Candies          4 31470.6  7867.7 780.1762 < 2e-16 ***
## Judges          10   224.9    22.5   2.2304 0.02089 *  
## Candies:Judges  40   707.7    17.7   1.7545 0.01158 *  
## Residuals      110  1109.3    10.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC  2                Df Sum Sq Mean Sq F value    Pr(>F)    
## Candies          4 1573.8  393.46 33.1604 < 2.2e-16 ***
## Judges          10  278.3   27.83  2.3455 0.0150274 *  
## Candies:Judges  40 1053.3   26.33  2.2193 0.0005888 ***
## Residuals      110 1305.2   11.87                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC  3                Df  Sum Sq Mean Sq F value    Pr(>F)    
## Candies          4  307.12  76.780   7.646 1.790e-05 ***
## Judges          10 1006.62 100.662  10.024 8.574e-12 ***
## Candies:Judges  40  484.02  12.101   1.205    0.2229    
## Residuals      110 1104.61  10.042                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC  4                Df Sum Sq Mean Sq F value Pr(>F)
## Candies          4  29.62  7.4041  1.0131 0.4039
## Judges          10 108.27 10.8275  1.4816 0.1557
## Candies:Judges  40 357.43  8.9357  1.2227 0.2062
## Residuals      110 803.89  7.3081               
## PC  5                Df Sum Sq Mean Sq F value  Pr(>F)  
## Candies          4  13.23  3.3069  0.7059 0.58958  
## Judges          10  62.48  6.2485  1.3338 0.22151  
## Candies:Judges  40 296.16  7.4040  1.5805 0.03254 *
## Residuals      110 515.32  4.6847                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC  6                Df Sum Sq Mean Sq F value   Pr(>F)   
## Candies          4  10.99  2.7468  0.5749 0.681450   
## Judges          10 133.09 13.3088  2.7852 0.004158 **
## Candies:Judges  40 176.20  4.4050  0.9219 0.605543   
## Residuals      110 525.62  4.7784                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC  7                Df Sum Sq Mean Sq F value    Pr(>F)    
## Candies          4   7.00  1.7500  0.7524 0.5585229    
## Judges          10 102.49 10.2492  4.4064  3.29e-05 ***
## Candies:Judges  40 208.38  5.2094  2.2397 0.0005144 ***
## Residuals      110 255.86  2.3260                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC  8                Df  Sum Sq Mean Sq F value Pr(>F)
## Candies          4   3.594  0.8986  0.3141 0.8680
## Judges          10  32.565  3.2565  1.1383 0.3407
## Candies:Judges  40 135.848  3.3962  1.1871 0.2408
## Residuals      110 314.693  2.8608
# p-values 
pander("p-values")

p-values

pander(pvalues)
Candies Judges Candies:Judges
1.442e-32 0.02089 0.01158
1.495e-07 0.01503 0.0005888
0.0004731 8.574e-12 0.2229
0.5149 0.1557 0.2062
0.7742 0.2215 0.03254
0.6484 0.004158 0.6055
0.8521 3.29e-05 0.0005144
0.8989 0.3407 0.2408
pander("pvalues <= 0.05")

pvalues <= 0.05

pvalues <= 0.05
##      Candies Judges Candies:Judges
## [1,]    TRUE   TRUE           TRUE
## [2,]    TRUE   TRUE           TRUE
## [3,]    TRUE   TRUE          FALSE
## [4,]   FALSE  FALSE          FALSE
## [5,]   FALSE  FALSE           TRUE
## [6,]   FALSE   TRUE          FALSE
## [7,]   FALSE   TRUE           TRUE
## [8,]   FALSE  FALSE          FALSE
pander("Bonferroni corrected p-values")

Bonferroni corrected p-values

pval_corrected <- t(round(pvalues*(nPC*3),4))
pval_corrected[pval_corrected>1]=1
pander(pval_corrected)
Candies 0 0 0.0114 1 1 1 1 1
Judges 0.5014 0.3607 0 1 1 0.0998 8e-04 1
Candies:Judges 0.278 0.0141 1 1 0.781 1 0.0123 1
# bonferroni corrected
pander("corrected pvalues <= 0.05")

corrected pvalues <= 0.05

pvalues <= 0.05/(nPC*3)
##      Candies Judges Candies:Judges
## [1,]    TRUE  FALSE          FALSE
## [2,]    TRUE  FALSE           TRUE
## [3,]    TRUE   TRUE          FALSE
## [4,]   FALSE  FALSE          FALSE
## [5,]   FALSE  FALSE          FALSE
## [6,]   FALSE  FALSE          FALSE
## [7,]   FALSE   TRUE           TRUE
## [8,]   FALSE  FALSE          FALSE
# write.csv(round(pvalues,4), file = file.path(out_path, "pvaluesBC.csv"))

6 GLM decomposition with fixed effects

6.1 Model definition

directoryInput <- "LIB"

# R scripts from M. Thiel
directory <-file.path(paste0(directoryInput,"/CodeMThiel"))
source(file.path(directory, "permutationTest.R"))
source(file.path(directory, "matrixDecomposition.R"))

# formula
formula <- as.formula(outcomes ~  Candies +  Judges +  Candies*Judges) 

6.2 GLM decomposition

# Decomposition of Y into the effect matrices by GLM

resGLM <- matrixDecomposition(formula,outcomes,design)
modelMatrix <- sapply(resGLM$modelMatrixByEffect, function(x) x)

nparam <- sum(sapply(modelMatrix, function(x) dim(x)[2]))-1

# Building of the list of pure effects matrices

EffectMatGLM <- resGLM$effectMatrices[-1] # minus intercept

res <- vector(mode = "list")
res[[1]]  <- resGLM$residuals

EffectMatGLM <- c(EffectMatGLM, residuals=res) # plus residuals

pander("names(EffectMatGLM)")

names(EffectMatGLM)

names(EffectMatGLM)
## [1] "Candies"        "Judges"         "Candies:Judges" "residuals"
ModelTerms <- names(EffectMatGLM)

6.3 scores/loadings generic graphs

listgraphs <- list()
varASCA <- list()

for(i in 1:length(ModelTerms)) {
  
  # PCA on the pure effect matrices
  ascaSVD = SVDforPCA(EffectMatGLM[[i]])
  ascaSVD$scores=round(ascaSVD$scores,5)
  varASCA[[i]] <- ascaSVD$var
  
  # Scores plot 
  a = DrawScores(ascaSVD, type.obj = "PCA", drawNames = TRUE,
                 createWindow = F, 
                 main = paste0(ModelTerms[i],"score plot - ASCA-GLM"), 
                 color = as.factor(Judges), 
                 pch = as.factor(Candies), axes = c(1, 2),size=2.5)

  # Loadings plot 
  b = ScatterPlot(ascaSVD$loadings[,1], 
                         ascaSVD$loadings[,2], 
                         points_labs = rownames(ascaSVD$loadings),
                          cex.lab = 3)+
  geom_vline(xintercept = 0, lwd=0.1) + 
  geom_hline(yintercept = 0, lwd=0.1)+
  labs(title = paste0(ModelTerms[i],"- loadings plot - ASCA-GLM"),
       x=paste0("PC1 (", round(ascaSVD$var[1],2), "%)"), 
       y=paste0("PC2 (", round(ascaSVD$var[2],2), "%)"))
  
  listgraphs[[paste0(ModelTerms[i],"_scores")]] <- a
  listgraphs[[paste0(ModelTerms[i],"_loadings")]] <- b
}

listgraphs
## $Candies_scores

## 
## $Candies_loadings

## 
## $Judges_scores

## 
## $Judges_loadings

## 
## $`Candies:Judges_scores`

## 
## $`Candies:Judges_loadings`

## 
## $residuals_scores

## 
## $residuals_loadings

7 Parallel mixed modelling

###################################################
# Parallel mixed modelling
###################################################

### Coding the interaction -------------------
CandiesJudges <- model.matrix(~ Candies:Judges-1, 
                                 data = design)


# number of parameters for the interaction
num <- 1:dim(CandiesJudges)[2]


for (i in 1:dim(CandiesJudges)[2]){
 CandiesJudges[which(CandiesJudges[,i]!=0),i] = i
}  

CandiesJudges <- rowSums(CandiesJudges)
CandiesJudges <- as.factor(CandiesJudges)

# New design with the coded interaction
designInter <- cbind(design,CandiesJudges)


### parallel LMM 
#######################

form <- "~ Candies +  (1|Judges) +  (1|CandiesJudges)"

REML <- TRUE

# run parlmer -------------------

# Modification of the parlmer function (=> into parlmer_interaction) 
# to obtain 2 by 2 orthogonal effect matrices.
res.parlmer <- parlmer_interaction(design=designInter, outcomes, form, REML)


MM_full <- res.parlmer

pander("summary for the first response")

summary for the first response

summary(res.parlmer$merMod_obj$PC1)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 876.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7692 -0.5925 -0.0261  0.4895  4.5514 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  CandiesJudges (Intercept)  1.1056  1.0515  
##  Judges        (Intercept)  0.7998  0.8943  
##  Residual                  10.4944  3.2395  
## Number of obs: 165, groups:  CandiesJudges, 55; Judges, 11
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)  5.154e-15  3.692e-01    0.00
## Candies1    -1.684e+01  5.958e-01  -28.27
## Candies2     1.167e+01  5.958e-01   19.59
## Candies3     1.036e+01  5.958e-01   17.40
## Candies4     1.177e+01  5.958e-01   19.77
## 
## Correlation of Fixed Effects:
##          (Intr) Cands1 Cands2 Cands3
## Candies1  0.000                     
## Candies2  0.000 -0.179              
## Candies3  0.000 -0.179 -0.179       
## Candies4  0.000 -0.179 -0.179 -0.179
# save the results -------------------

RanModMatlist <- res.parlmer$RanModMatlist

# head(RanModMatlist$`1 | CandiesJudges`)
# head(RanModMatlist$`1 | Judges`)

FixedModMatlist <- res.parlmer$FixedModMatlist
# head(FixedModMatlist$Candies)
# head(FixedModMatlist$`(Intercept)`)

# Residuals sd error
Res_std_error_PC <- sapply(MM_full$merMod_obj, sigma)

# Residuals
Residuals_PC <- sapply(MM_full$merMod_obj, residuals)


### ranef_PC: Extract the modes of the random effects
ranef_PC <- sapply(MM_full$merMod_obj, function(x) unlist(ranef(x, condVar=FALSE)))

# recover accurate rownames and colnames of ranef_PC
list_rownam <- lapply(ranef(MM_full$merMod_obj$PC1, condVar=FALSE), rownames)

colnam <- paste0(names(ranef(MM_full$merMod_obj$PC1, condVar=FALSE)), 
                 lapply(ranef(MM_full$merMod_obj$PC1, condVar=FALSE), rownames))

colnam <- rownames(ranef_PC) #(?)

rownames(ranef_PC) <- colnam
rownames(ranef_PC) <- gsub("..(Intercept))", "", rownames(ranef_PC))

### fixef: Extract fixed-effects estimates
fixef_PC <- sapply(MM_full$merMod_obj, fixef)

### all fixed estimates and random predictions
cof_PC <- rbind(fixef_PC, ranef_PC)

### Extract Variance and Correlation Components
varcor_random_full <- sapply(MM_full$merMod_obj, VarCorr) # Std.Dev.
# varcor_random_full[,1]
# dat <- as.numeric(rbind(varcor_random_full))
# varcor_random_full_mat <- matrix(dat, nrow=2, 
#                                  dimnames = dimnames(varcor_random_full))
# 

### varcor_fixed
varcor_fixed_full <- sapply(MM_full$merMod_obj, function(x) sqrt(diag(vcov(x)))) # Std.Dev.

rownames(varcor_fixed_full) <- rownames(vcov(MM_full$merMod_obj[[1]]))


##### Names of the effects
fixNames <- MM_full$fixNames # names of fixed effects

ranNames <- MM_full$ranNames # names of random effects

8 Normality of the model residuals ?

library("car")
# pdf(file.path(out_path, "SDA_residuals.pdf"), height = 8, width = 8)
par(mfrow=c(4,4), mar=c(4,4,2,2))
for (i in 1:ncol(Residuals_PC)){
  qqPlot(Residuals_PC[,i], main = paste0("qqplot - PC",i), 
         ylab="Residuals quantiles")
  d = density(Residuals_PC[,i])
  hist(Residuals_PC[,i], main = paste0("Histogram - PC",i), xlab="")
}
# dev.off()

9 Effect matrix computation

###################################################
# Effect matrix
###################################################
## Computation

# fixed effects + intercept -----------------------------
dim1FixedModMad <- sapply(FixedModMatlist, function(x) dim(x)[2])
names_FixedEffects <- names(FixedModMatlist)
shortFixedNames <- gsub("[^A-z]", "", names_FixedEffects)

Xmat <- do.call(cbind, FixedModMatlist)

# colnames(Xmat) %in% rownames(fixef_PC)
Xmat <- Xmat[,rownames(fixef_PC)] # reorder colnames of Xmat

index <- cumsum(dim1FixedModMad)
k <- 1
Mfix <- vector("list", length=length(shortFixedNames))
Mfix_PC <- vector("list", length=length(shortFixedNames))
names(Mfix) <- names(Mfix_PC) <- shortFixedNames
for (i in 1:length(shortFixedNames)){
  XMfix = Xmat
  XMfix[,-c(k:index[i])] = 0
  Mfix_PC[[i]] = XMfix %*% fixef_PC
  Mfix[[i]] <- Mfix_PC[[i]]%*%t(spectra_PCA_loadings)
  k <-  index[i] + 1
}

M0 <- Mfix[[1]]

Mfix <- Mfix[-1]
Mfix_PC <- Mfix_PC[-1]


# random effects -----------------------------
dim1RandModMad <- sapply(RanModMatlist, function(x) dim(x)[2])
names_randomEffects <- names(RanModMatlist)
shortRandNames <- gsub("[^A-z]", "", names_randomEffects)
# rbind(rownames(RanModMatlist$`1 | Judges`),rownames(RanModMatlist$`1 | CandiesJudges`)) # ok
Zmat <- do.call(cbind, RanModMatlist)
colnames(Zmat) <- paste0(rep(shortRandNames, dim1RandModMad),colnames(Zmat))

# colnames(Zmat) %in% rownames(ranef_PC)

index <- cumsum(dim1RandModMad)
k <- 1
Mrand <- vector("list", length=length(ranNames))
Mrand_PC <- vector("list", length=length(ranNames))
names(Mrand) <- names(Mrand_PC) <- ranNames
for (i in 1:length(shortRandNames)){
  XMrand = Zmat
  XMrand[,-c(k:index[i])] = 0
  Mrand_PC[[i]] = XMrand %*% ranef_PC
  Mrand[[i]] <- Mrand_PC[[i]]%*%t(spectra_PCA_loadings)
  k <-  index[i] + 1
}

# Residuals  -----------------------------
Mres_PC <- Residuals_PC
Mres <- Mres_PC%*%t(spectra_PCA_loadings)
Mres_full <- Mres

10 Are all the estimated effect matrices orthogonal ?

cat("t(Mrand$CandiesJudges)%*%Mrand$Judges")
## t(Mrand$CandiesJudges)%*%Mrand$Judges
pander(t(Mrand$CandiesJudges)%*%Mrand$Judges)
Table continues below
  Transp Acid Sweet Raspb Sugar
Transp 1.582e-15 -1.11e-16 -9.992e-15 -1.232e-14 1.277e-15
Acid 5.794e-15 9.176e-14 -1.948e-14 2.168e-14 9.284e-15
Sweet 4.108e-15 4.441e-16 -2.043e-14 -3.109e-15 4.885e-15
Raspb 3.886e-16 8.438e-15 1.288e-14 1.621e-14 -3.886e-15
Sugar -1.429e-15 4.108e-15 9.992e-15 -1.01e-14 -8.826e-15
Bites 8.743e-16 -7.022e-15 -3.22e-15 -4.302e-16 2.567e-16
Hard -3.102e-15 9.714e-15 4.774e-15 1.221e-15 -6.106e-16
Elastic 1.443e-15 -7.327e-15 7.438e-15 1.332e-15 4.718e-16
Sticky 1.36e-15 -4.441e-16 -3.109e-15 -2.887e-15 -3.331e-16
  Bites Hard Elastic Sticky
Transp -5.551e-16 -5.551e-15 -1.998e-15 -3.331e-16
Acid -4.385e-15 -1.638e-15 1.513e-15 -7.216e-16
Sweet 1.954e-14 -4.441e-16 -2.22e-16 -1.288e-14
Raspb 2.22e-15 -4.885e-15 -1.665e-15 7.55e-15
Sugar 3.331e-16 -3.22e-15 -1.221e-15 2.109e-15
Bites 2.914e-15 3.025e-15 -1.735e-16 2.359e-15
Hard -5.607e-15 -5.551e-17 1.693e-15 -2.387e-15
Elastic -4.774e-15 4.108e-15 -1.055e-15 -2.554e-15
Sticky 6.883e-15 6.661e-16 2.554e-15 5.107e-15
# unique(Ma[,1])
# unique(Mca[,1])
# 
# colSums(Ma)
# colSums(Mca)


cat("t(Mfix$Candies)%*%(Mrand$CandiesJudges)")
## t(Mfix$Candies)%*%(Mrand$CandiesJudges)
pander(t(Mfix$Candies)%*%(Mrand$CandiesJudges))
Table continues below
  Transp Acid Sweet Raspb Sugar
Transp 8.846e-13 1.088e-13 -1.3e-12 -1.03e-12 -7.994e-13
Acid -8.651e-13 -6.262e-14 1.116e-12 8.917e-13 7.656e-13
Sweet -2.269e-13 -5.052e-14 3.784e-13 3.046e-13 1.852e-13
Raspb -2.891e-13 -2.909e-14 4.032e-13 3.411e-13 2.42e-13
Sugar -9.273e-13 -1.292e-13 1.386e-12 1.073e-12 8.562e-13
Bites 8.278e-13 1.181e-13 -1.201e-12 -9.45e-13 -7.301e-13
Hard 8.278e-13 9.992e-14 -1.187e-12 -9.486e-13 -7.354e-13
Elastic 8.153e-13 1.137e-13 -1.208e-12 -9.486e-13 -7.407e-13
Sticky 7.327e-13 1.035e-13 -1.03e-12 -8.562e-13 -6.537e-13
  Bites Hard Elastic Sticky
Transp 8.786e-13 8.766e-13 7.514e-13 7.496e-13
Acid -8.484e-13 -8.695e-13 -7.621e-13 -7.532e-13
Sweet -2.42e-13 -2.338e-13 -1.779e-13 -1.767e-13
Raspb -2.874e-13 -2.9e-13 -2.3e-13 -2.549e-13
Sugar -9.213e-13 -9.335e-13 -7.834e-13 -7.816e-13
Bites 8.298e-13 8.251e-13 6.883e-13 6.892e-13
Hard 8.096e-13 8.251e-13 6.99e-13 7.141e-13
Elastic 8.107e-13 8.162e-13 6.883e-13 7.141e-13
Sticky 7.427e-13 7.443e-13 6.102e-13 6.217e-13

11 PCA on the effect matrices

11.1 PCA on the pure effect matrices

###################################################
## PCA on the effect matrices
###################################################


## Fixed effects -------------------------------
loadingsFix <- vector(mode = "list", length = length(Mfix))
for (i in 1:length(Mfix)){
  print(names(Mfix[i]))
  dimnames(Mfix[[i]]) <- dimnames(raw_outcomes)
  res_pca <- SVDforPCA(Mfix[[i]])
  res_pca$scores <- round(res_pca$scores, 6)
  
  barplot(res_pca$var[1:10], main  = "scree plot")
   
  print(DrawScores(res_pca, size=3, color = design[,names(Mfix[i])],
                         main = paste("Scores plot"),
                   axes=c(1,2), drawNames=FALSE))

  print(DrawScores(res_pca, size=3, color = design[,names(Mfix[i])],
             main = paste("Scores plot"), axes=c(3,4), 
        drawNames=FALSE))
  
  load <- DrawLoadings(res_pca, type.obj = "PCA", createWindow = F,
                       main = paste("Loadings plot"),
                       axes = c(1:4), loadingstype = "s",
                       num.stacked = 4, nxaxis = 9,
                       ang = "0",
                       xaxis_type = "character")
  
    load2 <- ScatterPlot(res_pca$loadings[,1], 
                         res_pca$loadings[,2], 
                         points_labs = rownames(res_pca$loadings),
                          cex.lab = 5)
    
    load3 <- ScatterPlot(res_pca$loadings[,3], 
                         res_pca$loadings[,4], 
                         points_labs = rownames(res_pca$loadings),
                          cex.lab = 5)
    print(load[[1]])
    print(load2)
    print(load3)
    loadingsFix[[i]] <-  load
}
## [1] "Candies"

## Random effects -------------------------------
loadingsRand <- vector(mode = "list", length = length(Mrand))
for (i in 1:length(Mrand)){
  print(names(Mrand[i]))
  dimnames(Mrand[[i]]) <- dimnames(raw_outcomes)
  res_pca <- SVDforPCA(Mrand[[i]])
  
  barplot(res_pca$var[1:10], main  = "scree plot")

  res_pca$scores <- round(res_pca$scores, 6)
  print(DrawScores(res_pca, size=3, color = design$Judges,
                   pch=design$Candies,
                   main = paste("Scores plot"),
                   axes=c(1,2), drawNames=FALSE))

  print(DrawScores(res_pca, size=3, color = design$Judges,
                   pch=design$Candies,
                   main = paste("Scores plot"), axes=c(3,4), 
                   drawNames=FALSE))
  
 
  load <- DrawLoadings(res_pca, axes=c(1:4), type.obj = "PCA",
                       xaxis_type = "character",
                       loadingstype = "s",  xlab = "ppm", nxaxis = 9,
                       main = paste("Loadings plot:"))
    
  load2 <- ScatterPlot(res_pca$loadings[,1], 
                       res_pca$loadings[,2], 
                       points_labs = rownames(res_pca$loadings),
                          cex.lab = 5)

   
    load3 <- ScatterPlot(res_pca$loadings[,3], 
                         res_pca$loadings[,4], 
                         points_labs = rownames(res_pca$loadings),
                          cex.lab = 5)
    print(load[[1]])
    print(load2)
    print(load3)
  loadingsRand[[i]] <-  load
}
## [1] "CandiesJudges"

## [1] "Judges"

11.1.1 Scree plots and pure loadings

Note that only the loadings are backtransformed

# ========================
# Pure loadings ==========
# ========================

# Candies
#####################
res_pca <- SVDforPCA(Mfix_PC$Candies)
colnames(Mfix_PC$Candies) <- paste0("Var",colnames(Mfix_PC$Candies))
df <- data.frame(PC = as.character(1:8), var = res_pca$var)

screeplotCandies <- ggplot(df, aes(y=0,yend=var,x = PC,
                          xend=PC))+ geom_segment()+
                   labs(title= "Scree plot",
                        x = "PC", y="% var") + theme_classic()
           
backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
loadCandies <- ScatterPlot(backtransf_load[,1], 
                         backtransf_load[,2], 
                         points_labs = rownames(backtransf_load),
                          cex.lab = 3)+
  geom_vline(xintercept = 0, lwd=0.1) + 
  geom_hline(yintercept = 0, lwd=0.1)+
  coord_cartesian(xlim = c(-0.5, 0.5), ylim = c(-1, 1))+
  labs(title = "PCA loadings", x=paste0("PC1 (", round(res_pca$var[1],2), "%)"),
       y=paste0("PC2 (", round(res_pca$var[2],2), "%)"))



## Random effects -------------------------------

# Judges ++++++++
res_pca <- SVDforPCA(Mrand_PC$Judges)
df <- data.frame(PC = as.character(1:length(res_pca$var)),
                 var = res_pca$var)

screeplotJudges <- ggplot(df, aes(y=0,yend=var,x=PC,
                          xend=PC))+ geom_segment()+
              labs(title= "Scree plot",
                        x = "PC", y="% var")+theme_classic()

  
# screeplot 
backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
loadJudges <- ScatterPlot(backtransf_load[,1], 
                         backtransf_load[,2], 
                         points_labs = rownames(backtransf_load),
                          cex.lab = 3)+
  geom_vline(xintercept = 0, lwd=0.1) + 
  geom_hline(yintercept = 0, lwd=0.1)+
  coord_cartesian(xlim = c(-0.8, 0.8), ylim = c(-0.6, 0.6))+
  labs(title = "PCA loadings", x=paste0("PC1 (", round(res_pca$var[1],2), "%)"), y=paste0("PC2 (", round(res_pca$var[2],2), "%)"))


# CandiesJudges ++++++++
res_pca <- SVDforPCA(Mrand_PC$CandiesJudges)
df <- data.frame(PC = as.character(1:length(res_pca$var)),
                 var = res_pca$var)

screeplotCandiesJudges <- ggplot(df, aes(y=0,yend=var,x=PC,
                          xend=PC))+ geom_segment()+
          labs(title= "Scree plot",
                        x = "PC", y="% var")+theme_classic()

# screeplot 
backtransf_load <- -1*t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
loadCandiesJudges <- ScatterPlot(backtransf_load[,1], 
                         backtransf_load[,2], 
                         points_labs = rownames(backtransf_load),
                          cex.lab = 3) +
  geom_vline(xintercept = 0, lwd=0.1) + 
  geom_hline(yintercept = 0, lwd=0.1)+
  coord_cartesian(xlim = c(-0.8, 0.8), ylim = c(-0.6, 0.6)) +
  labs(title = "PCA loadings", x=paste0("PC1 (", round(res_pca$var[1],2), "%)"), y=paste0("PC2 (", round(res_pca$var[2],2), "%)"))

11.2 PCA on the augmented effect matrices (ASCA-E)

11.2.1 addSegments fun

addSegments <- function(group, data, pch=16, main = NULL, 
                        col = rainbow(n = length(unique(group))),
                        ...) {
  
  group <- as.factor(group)
  
  # col <- rainbow(n = length(unique(group)))
  
  plot(data$x, data$y, col=col[group], pch=pch, main=main, ...)
  
  group <- as.factor(group)
  xcent <- tapply(data[,1], group, FUN=mean)
  ycent <- tapply(data[,2], group, FUN=mean)
  centers <- data.frame(xcent=xcent, ycent=ycent)
  
  mapply(FUN = points, x = centers$xcent, y = centers$ycent, 
         col = col, MoreArgs = list(pch=20, cex=0.7))
  
  submatrices <- split(x=data, f=group)
  for (i in 1:nlevels(group)){
    mapply(FUN = segments, x1 = submatrices[[i]]$x, 
           y1 = submatrices[[i]]$y, col = col[i],
           MoreArgs = list(x0 = centers$xcent[i], 
                           y0 = centers$ycent[i]))
  }
  
  
}

11.2.2 Augmented (ANOVA) scores plots

# ========================
# augmented scores =======
# ========================
# * Candies: Mc + Mca
# * Assessor: Ma + E
# * Candies:Assessor : Mca + E

a <- nlevels(design$Candies)
b <- nlevels(design$Judges)
nn <- table(design$Judges,design$Candies)[1,1]

# Candies
##################
# pander("Candies")
ascaSVD = SVDforPCA(Mfix_PC$Candies)
Fstat <- qf(.95, df1=(a-1), df2=((a-1)*(b-1))) 
coef <- sqrt(Fstat/(b-1))
ascaSVD$scores[,1:2]=(Mfix_PC$Candies + 
                        ((Mrand_PC$CandiesJudges)*coef)) %*%
                        ascaSVD$loadings[,1:2]


Candies_scores <- DrawScores(ascaSVD, 
                main = "PCA Scores",
           color = Candies, pch = Candies, 
           drawNames = FALSE, drawPolygon  = FALSE, drawEllipses = TRUE,
           noLegend = FALSE, size=2) 


# Judges
##################
# pander("Judges")
ascaSVD = SVDforPCA(Mrand_PC$Judges)
df1 <- (b-1)
df2 <- (a*b*(nn-1))
Fstat <- qf(.95, df1=df1, df2=df2) 
coef <- sqrt(Fstat*df1/df2)

ascaSVD$scores[,1:2]=(Mrand_PC$Judges+(Mres_PC*coef))%*%
            ascaSVD$loadings[,1:2]
Judges_scores <- DrawScores(ascaSVD, type.obj = "PCA", 
                drawNames = F, createWindow = F, 
           main = "PCA Scores", 
           color = Judges, 
           pch = Judges, axes = c(1, 2),size=2, drawPolygon = FALSE, 
           drawEllipses = TRUE)+ 
  theme(legend.text=element_text(size=10),
        legend.key.height=unit(0.7,"line"))




# CandiesJudges
##################
# pander("CandiesJudges")
ascaSVD = SVDforPCA(Mrand_PC$CandiesJudges)

df1 <- (a-1)*(b-1)
df2 <- (a*b*(nn-1))
Fstat <- qf(.95, df1=df1, df2=df2) 
coef <- sqrt(Fstat*df1/df2)

ascaSVD$scores[,1:2]=-1*(Mrand_PC$CandiesJudges+(Mres_PC*coef))%*%
          ascaSVD$loadings[,1:2]



CJ <- as.factor(paste0(design$Candies, design$Judges))
index <- c(145,147,27, 109, 110)

data <- ascaSVD$scores[,1:2]
group <- as.factor(CJ)
xcent <- tapply(data[,1], group, FUN=mean)
ycent <- tapply(data[,2], group, FUN=mean)
xc <- yc <- c()
for (i in 1:length(group)){
  xc[i] <- xcent[group[i]]
  yc[i] <- ycent[group[i]]
}

df <- data.frame(data, xcent=xc, ycent = yc, CJ=CJ,
                 Judges =design$Judges, Candies=design$Candies)

# Eigenvalues
eig <- ascaSVD$eigval
# Variances in percentage
variance <- eig * 100/sum(eig)
Xax=1
Yax=2
XaxName <- paste0("PC", Xax, " (", round(variance[Xax], 2),"%)")
YaxName <- paste0("PC", Yax, " (", round(variance[Yax], 2), "%)")
xlab <- XaxName
ylab <- YaxName
Xlim <- c(min(ascaSVD$scores[, Xax]) * 1.4, max(ascaSVD$scores[, Xax]) * 1.4)
Ylim <- c(min(ascaSVD$scores[, Yax]) * 1.4, max(ascaSVD$scores[, Yax]) * 1.4)
main = "PCA Scores"

plots <- ggplot(df,aes(x=xc,y=yc))+
  geom_point(df, mapping=aes(x=PC1,y=PC2, shape=Candies, color=Judges),size=2) +
  geom_segment(aes(yend=PC2,xend=PC1,color=Judges,group=CJ)) + 
  ggplot2::xlim(Xlim) + ggplot2::ylim(Ylim)  +
            scale_shape_manual(values=seq(0,26), name = "Candies")

plots <- plots + ggplot2::labs(title = main, x = xlab, y = ylab) + 
  ggplot2::geom_vline(xintercept = 0,size = 0.1) + 
  ggplot2::geom_hline(yintercept = 0, size = 0.1) + ggplot2::theme_bw() +
  ggplot2::theme(panel.grid.major = 
                   ggplot2::element_line(color = "gray60", size = 0.2), 
                 panel.grid.minor = ggplot2::element_blank(),
                 panel.background = ggplot2::element_rect(fill = "gray98"))

plots <- plots + theme(legend.text=element_text(size=10),
                       legend.key.height=unit(0.6,"line")) +
  annotate("text", y = (ascaSVD$scores[index,2] +
                                1.7*c(-1 ,1,-1,-1, -1)),
                         x = ascaSVD$scores[index,1],
                         label = row.names(outcomes[index,1:2]))

CA_scores <- plots



mat <- data.frame(x=ascaSVD$scores[,1],y=ascaSVD$scores[,2])
group <- as.factor(paste0(Candies,Judges))
# table(group)

# pdf(file.path(fig_path,"SDA_ScoresCJEM_linkbyCJ.pdf"))
# addSegments(group=group, data=mat, main = "Scores plot C*J effect matrix \n color by CandiesJudges", ylab = "PC2", xlab = "PC1")
# abline(h=0, v=0, lty=2, lwd =1,col="black")
# dev.off()


# Residuals ================

res_pca <- SVDforPCA(Mres_PC)
# colnames(Mfix_PC$Candies) <- paste0("Var",colnames(Mfix_PC$Candies))
df <- data.frame(PC = as.character(1:8), var = res_pca$var)
# df <- df[1:7,]
screeplot_resid <- ggplot(df, aes(y=0,yend=var,x = PC,
                          xend=PC))+ geom_segment()+
                   labs(title= "Scree plot",
                        x = "PC", y="% var") + theme_classic()
           
backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))

loadings_resid <- ScatterPlot(backtransf_load[,1], 
                         backtransf_load[,2], 
                         points_labs = rownames(backtransf_load),
                          cex.lab = 3)+
  geom_vline(xintercept = 0, lwd=0.1) + 
  geom_hline(yintercept = 0, lwd=0.1)+
  labs(title = "PCA  loadings", 
       x=paste0("PC1 (", round(res_pca$var[1],2), "%)"), y=paste0("PC2 (", round(res_pca$var[2],2), "%)")) +
  coord_cartesian(xlim = c(-0.85, 0.85), ylim = c(-0.5, 0.5))

index <- c(10, 27)

df <- data.frame(PC = as.character(1:length(res_pca$var)),
                 var = res_pca$var)
df <- df[1:7,]

scores_resid <- DrawScores(res_pca, drawNames = FALSE, 
                           color = Candies, 
           pch = Judges, main ="PCA scores",
           size = 2) +
  coord_cartesian(xlim = c(-16, 16), ylim = c(-18, 18))+ 
  theme(legend.text=element_text(size=10), legend.key.height=unit(0.7,"line"))+
  annotate("text", y = (res_pca$scores[index,2] +
                                1.7*c(-1 , -1)),
                         x = res_pca$scores[index,1],
                         label = rownames(res_pca$scores[index,1:2]))

11.2.3 Augmented (ED) scores plots

pander("Effective dimensions")

Effective dimensions

# Effective Dimensions
source("ED/Sensory_Mixed.R")

EDc <- rep((nlevels(design$Candies)-1),ncol(outcomes))
EDj <- ED_sensory["ED_judges",]
EDcj <- ED_sensory["ED_cj",]
EDe <- n - colSums(ED_sensory)
  
# ========================
# augmented scores =======
# ========================
# * Candies: Mc + Mca
# * Assessor: Ma + E
# * Candies:Assessor : Mca + E


# Candies
##################
# pander("Candies")
ascaSVD = SVDforPCA(Mfix_PC$Candies)
df1 <- EDc
df2= EDcj
Fstat <- qf(.95, df1=df1 ,df2= df2)
coef <- sqrt((Fstat*df1)/df2)

mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
for (i in 1:ncol(outcomes)){
  mat[,i] <- Mrand_PC$CandiesJudges[,i]*coef[i]
}

ascaSVD$scores[,1:2]= (Mfix_PC$Candies + mat) %*%
                        ascaSVD$loadings[,1:2]


Candies_scores <- DrawScores(ascaSVD, 
                main = "PCA Scores",
           color = Candies, pch = Candies, 
           drawNames = FALSE, drawPolygon  = FALSE, drawEllipses = TRUE,
           noLegend = FALSE, size=2) 


# Judges
##################
# pander("Judges")
ascaSVD = SVDforPCA(Mrand_PC$Judges)
df1 <- EDj
df2 <- EDe
Fstat <- qf(.95, df1=df1, df2=df2) 
coef <- sqrt(Fstat*df1/df2)

mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
for (i in 1:ncol(outcomes)){
  mat[,i] <- Mres_PC[,i]*coef[i]
}

ascaSVD$scores[,1:2]=(Mrand_PC$Judges + mat)%*%
            ascaSVD$loadings[,1:2]
Judges_scores <- DrawScores(ascaSVD, type.obj = "PCA", 
                drawNames = F, createWindow = F, 
           main = "PCA Scores", 
           color = Judges, 
           pch = Judges, axes = c(1, 2),size=2, drawPolygon = FALSE, 
           drawEllipses = TRUE)+ 
  theme(legend.text=element_text(size=10),
        legend.key.height=unit(0.7,"line"))



# CandiesJudges
##################
# pander("CandiesJudges")
ascaSVD = SVDforPCA(Mrand_PC$CandiesJudges)

df1 <- EDcj
df2 <- EDe
Fstat <- qf(.95, df1=df1, df2=df2) 
coef <- sqrt(Fstat*df1/df2)

mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
for (i in 1:ncol(outcomes)){
  mat[,i] <- Mres_PC[,i]*coef[i]
}


ascaSVD$scores[,1:2]=-1*(Mrand_PC$CandiesJudges + mat)%*%
          ascaSVD$loadings[,1:2]



CJ <- as.factor(paste0(design$Candies, design$Judges))
index <- c(145,147,27, 109, 110)

data <- ascaSVD$scores[,1:2]
group <- as.factor(CJ)
xcent <- tapply(data[,1], group, FUN=mean)
ycent <- tapply(data[,2], group, FUN=mean)
xc <- yc <- c()
for (i in 1:length(group)){
  xc[i] <- xcent[group[i]]
  yc[i] <- ycent[group[i]]
}

df <- data.frame(data, xcent=xc, ycent = yc, CJ=CJ,
                 Judges =design$Judges, Candies=design$Candies)

# Eigenvalues
eig <- ascaSVD$eigval
# Variances in percentage
variance <- eig * 100/sum(eig)
Xax=1
Yax=2
XaxName <- paste0("PC", Xax, " (", round(variance[Xax], 2),"%)")
YaxName <- paste0("PC", Yax, " (", round(variance[Yax], 2), "%)")
xlab <- XaxName
ylab <- YaxName
Xlim <- c(min(ascaSVD$scores[, Xax]) * 1.4, max(ascaSVD$scores[, Xax]) * 1.4)
Ylim <- c(min(ascaSVD$scores[, Yax]) * 1.4, max(ascaSVD$scores[, Yax]) * 1.4)
main = "PCA Scores"

plots <- ggplot(df,aes(x=xc,y=yc))+
  geom_point(df, mapping=aes(x=PC1,y=PC2, shape=Candies, color=Judges),size=2) +
  geom_segment(aes(yend=PC2,xend=PC1,color=Judges,group=CJ)) + 
  ggplot2::xlim(Xlim) + ggplot2::ylim(Ylim)  +
            scale_shape_manual(values=seq(0,26), name = "Candies")

plots <- plots + ggplot2::labs(title = main, x = xlab, y = ylab) + 
  ggplot2::geom_vline(xintercept = 0,size = 0.1) + 
  ggplot2::geom_hline(yintercept = 0, size = 0.1) + ggplot2::theme_bw() +
  ggplot2::theme(panel.grid.major = 
                   ggplot2::element_line(color = "gray60", size = 0.2), 
                 panel.grid.minor = ggplot2::element_blank(),
                 panel.background = ggplot2::element_rect(fill = "gray98"))

plots <- plots + theme(legend.text=element_text(size=10),
                       legend.key.height=unit(0.6,"line"))

CA_scores <- plots

# CA_scores

mat <- data.frame(x=ascaSVD$scores[,1],y=ascaSVD$scores[,2])
group <- as.factor(paste0(Candies,Judges))
# table(group)

# pdf(file.path(fig_path,"SDA_ScoresCJEM_linkbyCJ.pdf"))
# addSegments(group=group, data=mat, main = "Scores plot C*J effect matrix \n color by CandiesJudges", ylab = "PC2", xlab = "PC1")
# abline(h=0, v=0, lty=2, lwd =1,col="black")
# dev.off()

# ========================
# Residuals ==============
# ========================


res_pca <- SVDforPCA(Mres_PC)
# colnames(Mfix_PC$Candies) <- paste0("Var",colnames(Mfix_PC$Candies))
df <- data.frame(PC = as.character(1:8), var = res_pca$var)
# df <- df[1:7,]
screeplot_resid <- ggplot(df, aes(y=0,yend=var,x = PC,
                          xend=PC))+ geom_segment()+
                   labs(title= "Scree plot",
                        x = "PC", y="% var") + theme_classic()
           
backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))

loadings_resid <- ScatterPlot(backtransf_load[,1], 
                         backtransf_load[,2], 
                         points_labs = rownames(backtransf_load),
                          cex.lab = 3)+
  geom_vline(xintercept = 0, lwd=0.1) + 
  geom_hline(yintercept = 0, lwd=0.1)+
  labs(title = "PCA  loadings", 
       x=paste0("PC1 (", round(res_pca$var[1],2), "%)"), y=paste0("PC2 (", round(res_pca$var[2],2), "%)")) +
  coord_cartesian(xlim = c(-0.85, 0.85), ylim = c(-0.5, 0.5))

index <- c(10, 27)

df <- data.frame(PC = as.character(1:length(res_pca$var)),
                 var = res_pca$var)
df <- df[1:7,]

scores_resid <- DrawScores(res_pca, drawNames = FALSE, 
                           color = Candies, 
           pch = Judges, main ="PCA scores",
           size = 2) +
  coord_cartesian(xlim = c(-19, 19), ylim = c(-18, 18))+ 
  theme(legend.text=element_text(size=10), legend.key.height=unit(0.7,"line"))+
  annotate("text", y = (res_pca$scores[index,2] +
                                1.7*c(-1 , -1)),
                         x = res_pca$scores[index,1],
                         label = rownames(res_pca$scores[index,1:2]))

11.2.4 Graphs

Scores_Loadings_EffectMat <- ggarrange(screeplotCandies, 
                                       Candies_scores,
                                       loadCandies, 
                                       screeplotJudges,
                                       Judges_scores,
                                       loadJudges, 
                                       screeplotCandiesJudges, 
                                       CA_scores,
                                       loadCandiesJudges,
          labels = c("A", "B", "C", "D", "E", "F", "G", "H", "I"),
          ncol = 3, nrow = 3, common.legend=TRUE, widths = c(0.8,2,2))

# Scores_Loadings_EffectMat

# ggexport(Scores_Loadings_EffectMat, filename = file.path(fig_path,"SDA_Scores_Loadings_EffectMat.pdf"),
#          height = 12, width = 12)


a <- grid.arrange(screeplotCandies, Candies_scores,
              loadCandies, nrow=1,widths=c(0.3, 1,  0.85),
             top=textGrob("Candy effect matrix",
                          gp=gpar(fontsize=20,font=2)))

b <-  grid.arrange(screeplotJudges, Judges_scores,
               loadJudges, nrow=1,widths=c(0.3, 1,  0.85),
             top=textGrob("Judge effect matrix",
                          gp=gpar(fontsize=20,font=2)))

c <- grid.arrange(screeplotCandiesJudges, 
              CA_scores, loadCandiesJudges,
              nrow=1,widths  = c(0.3, 1,  0.85),
             top=textGrob("C*J effect matrix",
                          gp=gpar(fontsize=20,font=2)))

d <- grid.arrange(screeplot_resid,
              scores_resid, loadings_resid,
              nrow=1,widths  = c(0.3, 1,  0.85),
             top=textGrob("Residual effect matrix",
                          gp=gpar(fontsize=20,font=2)))

# Thesis chapter output ----------


p <- ggarrange(a,b,c, nrow=3, labels = c("A", "B", "C"))

# ggsave(file.path(fig_path,"SDA_Scores_Loadings_EffectMat.pdf"),
#          plot = p, height = 15, width = 14,scale = 0.65)



# journal article output ----------


plots <- ggarrange(a,b,c,d,  nrow=4, labels = c("A", "B", "C", "D"))
# ggsave(file.path(fig_path,"SDA_Scores_Loadings_EffectMat.pdf"), 
#          plot = plots,  height = 20, width = 14, scale = 0.7)

11.2.5 Spotted outliers

# outliers detected
index <- which(design$Judges=="07" & design$Candies=="1")
tab_a <- raw_outcomes[index,]
rownames(tab_a) <- paste0("07_1_",1:3)

index <- which(design$Judges=="02" & design$Candies=="1")
tab_b <- raw_outcomes[index,]
rownames(tab_b) <- paste0("02_1_",1:3)

# ggexport(ggtexttable(tab_a, theme = ttheme("classic")), 
#          filename =
#            file.path(fig_path,"SDA_Residuals_outcomes_07_1.pdf"),
#          height = 2, width = 6)

# ggexport(ggtexttable(tab_b, theme = ttheme("classic")), 
#          filename =
#            file.path(fig_path,"SDA_Residuals_outcomes_02_1.pdf"),
#          height = 2, width = 6)

11.2.6 Interaction graphs

11.2.6.1 Pure scores

# 1. Pure scores
#############################

res_pca <- SVDforPCA(Mrand_PC$CandiesJudges)
id <- which(!duplicated(res_pca$scores[,1]))
df <- as.data.frame(cbind(Candies=Candies[id], 
                 Judges = Judges[id],
                 scores = res_pca$scores[id,]))


df$Candies <- as.factor(df$Candies)
df$Judges <- as.factor(df$Judges)

a <- ggplot(data=df, aes(x=Candies, y=PC1)) + 
  geom_point(aes(colour = Judges, shape = Judges), size=3) +  
  scale_shape_manual(values = c(0:10)) + 
  geom_hline(yintercept=0, linetype=3, color = "black")

b <- ggplot(data=df, aes(x=Judges, y=PC1)) + 
  geom_point(aes(colour = Candies, shape =Candies), size=3)+ 
  geom_hline(yintercept=0, linetype=3, color = "black")


c <- ggplot(data=df, aes(x=Candies, y=PC2)) + 
  geom_point(aes(colour = Judges, shape = Judges), size=3) +  
  scale_shape_manual(values = c(0:10))+ 
  geom_hline(yintercept=0, linetype=3, color = "black")

d <-ggplot(data=df, aes(x=Judges, y=PC2))+
  geom_point(aes(colour = Candies, shape =Candies), size=3)+ 
  geom_hline(yintercept=0, linetype=3, color = "black")
# ggexport(ggarrange(a, c, common.legend=TRUE, legend="right"), 
#          filename = file.path(fig_path,"SDA_Scores_CAinteraction.pdf"),
#          height = 4, width = 12)

# pdf(file  = file.path(fig_path,"SDA_Scores_CAinteraction.pdf"),
#     height = 4, width = 12)
grid.arrange(a, c, nrow=1,widths=c(1,  0.85),
             top=textGrob("PCA scores on augmented C*J effect matrix",gp=gpar(fontsize=20,font=2)))

# dev.off()

11.2.6.2 Augmented scores

# 2. Augmented scores
#############################

# a <- nlevels(design$Candies)
# b <- nlevels(design$Judges)
# nn <- table(design$Judges,design$Candies)[1,1]

EDc <- rep((nlevels(design$Candies)-1),ncol(outcomes))
EDj <- ED_sensory["ED_judges",]
EDcj <- ED_sensory["ED_cj",]
EDe <- n - colSums(ED_sensory)

pander("CandiesJudges")

CandiesJudges

res_pca = SVDforPCA(Mrand_PC$CandiesJudges)

# df1 <- (a-1)*(b-1)
# df2 <- (a*b*(nn-1))
# Fstat <- qf(.95, df1=df1, df2=df2) 
# coef <- sqrt(Fstat*df1/df2)
# 
# res_pca$scores[,1:2]=-1*(Mrand_PC$CandiesJudges+(Mres_PC*coef))%*%
#           res_pca$loadings[,1:2]

df1 <- EDcj
df2 <- EDe
Fstat <- qf(.95, df1=df1, df2=df2) 
coef <- sqrt(Fstat*df1/df2)

mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
for (i in 1:ncol(outcomes)){
  mat[,i] <- Mres_PC[,i]*coef[i]
}


res_pca$scores[,1:2]=-1*(Mrand_PC$CandiesJudges + mat)%*%
          res_pca$loadings[,1:2]


df <- as.data.frame(cbind(Candies=Candies, 
                 Judges = Judges,
                 scores = res_pca$scores))


df$Candies <- as.factor(df$Candies)
df$Judges <- as.factor(df$Judges)

a <- ggplot(data=df, aes(x=Candies, y=PC1)) + geom_point(aes(colour = Judges, shape = Judges), size=3) +  
  scale_shape_manual(values = c(0:10)) + 
  geom_hline(yintercept=0, linetype=3, color = "black")+
  theme(text=element_text(size=15), 
        legend.key.height  = unit("0.7", units = "line"))

b <- ggplot(data=df, aes(x=Judges, y=PC1)) + geom_point(aes(colour = Candies, shape =Candies), size=3)+ 
  geom_hline(yintercept=0, linetype=3, color = "black")+
  theme(text=element_text(size=15), 
        legend.key.height  = unit("0.7", units = "line"))
# interPC1 <- grid.arrange(a, b)

c <- ggplot(data=df, aes(x=Candies, y=PC2)) + geom_point(aes(colour = Judges, shape = Judges), size=3) +  
  scale_shape_manual(values = c(0:10))+ 
  geom_hline(yintercept=0, linetype=3, color = "black")+
  theme(text=element_text(size=15), 
        legend.key.height  = unit("0.7", units = "line"))

d <-ggplot(data=df, aes(x=Judges, y=PC2))+geom_point(aes(colour = Candies, shape =Candies), size=3)+ 
  geom_hline(yintercept=0, linetype=3, color = "black")+
  theme(text=element_text(size=15), 
        legend.key.height  = unit("0.7", units = "line"))
# interPC2 <- grid.arrange(c, d)


a

b

c

d

# ggexport(ggarrange(a, c, common.legend=TRUE, legend="right"), 
#          filename = file.path(fig_path,"SDA_Scores_CAinteraction_aug.pdf"),
#          height = 4, width = 12)


p <- grid.arrange(a, c, nrow=1,widths=c(1,  1),
             top=textGrob("PCA scores on augmented C*J effect matrix",
                          gp=gpar(fontsize=20,font=2)))

# ggsave(file.path(fig_path,"SDA_Scores_CAinteraction_aug.pdf"),
#        plot = p, height = 4, width = 12, scale=0.9)

12 Percentage of explained variance

12.1 Method from Nakagawa and Schielzeth (2012)

####################################################
# Percentage of explained variance
####################################################
## Method from Nakagawa and Schielzeth (2012)

# Random effects -----------------------------------
sigma2_res = Res_std_error_PC^2 # Residual
varcor_random_full <- as.data.frame(varcor_random_full)
Var_Mrand <- rbind(varcor_random_full, sigma2_res=sigma2_res) # only random effects
Var_Mrand <- data.matrix(Var_Mrand)



# fixed effect -----------------------------------
Var_Mfix <- c()
for (i in 1:length(fixNames)){
  #  variance of parameters values (population)
  Var_Mfix <- rbind(Var_Mfix, (apply(Mfix_PC[[i]], 2, var) *(n - 1) / n))
}


# all together -----------------------------------
rownames(Var_Mfix) <- fixNames
rownames(Var_Mrand) <- c(ranNames, "Residuals")
var_comp <- rbind(Var_Mfix, Var_Mrand)
var_comp <- var_comp[c(1,3,2,4),]

# log of var comp +++++
log_var_comp <- t(log1p(var_comp)) # log(x+1)
log_var_comp <- cbind(id=rownames(log_var_comp), log_var_comp)
log_var_comp <- as.data.frame(log_var_comp)
log_var_comp <- melt(log_var_comp, id=c("id"))
log_var_comp$value <- as.numeric(log_var_comp$value)

names(log_var_comp) <- c("PC", "Effect", "Variance")
  
sum_var_comp <- rowSums(var_comp)

# Percent var expl by each effect -----
var_components_abs <- sum_var_comp
var_components <- var_components_abs*100/sum(var_components_abs)
names(var_components) <- rownames(var_comp)


# pdf(file.path(fig.path,"variance_components.pdf"), height = 3, width = 3, pointsize = 12)
par(mar=c(0.5,2,3,0.5))
barplot(var_components, main="Variance components \n percentage",xaxt="n",las=2, col=c(darkblue,turquoise,violetred , limegreen, gray67), border = NA,
    legend = names(var_components), args.legend = list(x="topright", 
            inset=c(0,0),box.lty=0,cex = 1, y.intersp = 0.8))
# dev.off()


table_var <- rbind(var_components_abs, var_components)
rownames(table_var) <- c("Sum variance for all responses", "percentage of variation")
colnames(table_var) <- rownames(var_comp)
pander(table_var)
  Candies Judges CandiesJudges Residuals
Sum variance for all responses 202.5 9.442 5.644 53.24
percentage of variation 74.77 3.486 2.084 19.66
# write.csv(table_var, file =  file.path(fig_path,"SDA_table_var.csv"))

12.2 Graph

log_var_comp$PC <- sub("Var", "", log_var_comp$PC)

p <- ggplot(log_var_comp, aes(Effect, Variance, group = PC))+
  ggtitle("SDA - Log of variance components")

p <- p + geom_point(aes(colour = PC))+ 
  geom_line(aes(colour = PC, linetype=PC),size=0.5)+
  theme_classic()+
  theme(legend.key.width = unit(0.8,"cm"),
        plot.title = element_text(size = 17),
        axis.title.y = element_text(size = 17),
        axis.title.x = element_text(size = 17),
        axis.text = element_text(size = 13),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 17)) + 
  ylab(label = "log(variance)")


# Thesis chapter output
# ggexport(p, filename = file.path(fig_path,"SDA_variance_components.pdf"),
#          height = 5, width = 6.5, pointsize=20)

tab <- data.frame(Effect= names(var_components), 
                  pcvar = round(var_components,2))
var_comp.table <- ggtexttable(tab, cols = c("Effect", "Global var (%)"), 
                                 rows = NULL, 
                                 theme = ttheme("classic",base_size = 10))

p <- p + annotation_custom(ggplotGrob(var_comp.table),
                              xmin = 6, ymin = 3,
                              xmax = 0) + 
  theme(legend.title = element_blank())

# journal article output
# ggexport(p, filename = file.path(fig_path,"SDA_variance_components.pdf"),
#          height = 5, width = 6.5, pointsize=19)
p

13 Bootstrap

13.1 Set up

# set up of the bootstrap
set.seed(2018)
nsim = 2000 # number of simulations

# name of the output 
name_RData <- "bootstrap_Sensory_Data.RData"

# Set up -------------------
# formulas without the effet to test
null_formulas <- list(Candies = "~ (1|Judges) +  (1|CandiesJudges)",
                       Judges = "~ Candies +  (1|CandiesJudges)",
                       CandiesJudges = "~Candies +  (1|Judges)")

null_effects <- names(null_formulas)

REML <- c(FALSE, TRUE, TRUE)
names(REML) <- null_effects

13.2 log-likelihood Ratio computation

13.2.1 True log-likelihood Ratio statistics

##################################################
# True log-likelihood Ratio statistics      #
##################################################


# full model: MM_full
######################
# REML +++++
loglik_PC_full_REML <- sapply(MM_full$merMod_obj, logLik, REML=T)

# ML +++++
loglik_PC_full_ML <- sapply(MM_full$merMod_obj, logLik, REML=F)

loglik_PC_full <- matrix(NA, ncol = nPC,
                         nrow=length(null_formulas), byrow = TRUE)
for (i in 1:length(REML)){
  if (REML[i]==TRUE){
    loglik_PC_full[i,] <- loglik_PC_full_REML
  }else {loglik_PC_full[i,] <- loglik_PC_full_ML}
}




### Restricted models  
######################
res.parlmer_NULL <- vector("list", length = length(null_formulas))
names(res.parlmer_NULL) <- names(REML) <- names(null_formulas)

for (i in 1:length(null_formulas)) {
  # run parlmer
  res.parlmer_NULL[[i]] <- parlmer(designInter, 
                                   outcomes, null_formulas[[i]], REML=REML[i])
}

# Save the results ------------
MM_PC_null <- lapply(res.parlmer_NULL, function(x) x[["merMod_obj"]])

#### Randvarnames, Fixvarnames
ranNames <- sapply(res.parlmer_NULL, function(x) x[["ranNames"]])
Fixvarnames <- sapply(res.parlmer_NULL, function(x) x[["fixNames"]])


varcor_random <- vector(mode = "list", length = length(null_formulas))
fixef_PC <- vector(mode = "list", length = length(null_formulas))
modmat_fixed <- vector(mode = "list", length = length(null_formulas))
M0 <- vector(mode = "list", length = length(null_formulas))

for (i in 1:length(null_formulas)){
  varcor_random[[i]] <- sapply(MM_PC_null[[i]], function(x)
    as.data.frame(VarCorr(x))$vcov)
  rownames(varcor_random[[i]]) <- c(names(VarCorr(MM_PC_null[[i]][[1]])),
                                    "Residual")

  fixef_PC[[i]] <- sapply(MM_PC_null[[i]], fixef)

  modmat_fixed[[i]] <- model.matrix(MM_PC_null[[i]][[1]], type = "fixed")

  # intercept
  XM0 <- modmat_fixed[[i]]
  XM0[,-1] <- 0
  M0_PC <- XM0%*%fixef_PC[[i]]
  M0[[i]] <- M0_PC%*%t(spectra_PCA_loadings)
}

# Effect matrices computation ------------
index <- vector("list", length=length(null_formulas))
fixNames_int <- lapply(Fixvarnames, function(x) c("(Intercept)", x))

colnam <- sapply(modmat_fixed,  colnames)

index <- vector("list", length=length(Fixvarnames))

for (i in 1:length(null_formulas)){
  id <- c()
  for (k in 1:length(fixNames_int[[i]])){
    id <- c(id, grep(fixNames_int[[i]][[k]], colnam[[i]]))
  }
  index[[i]] <- id
}


Mfix <- Mfix_PC <- vector("list", length=length(null_formulas))
names(Mfix) <- fixNames_int

for (i in 1:length(null_formulas)){
  XMfix = modmat_fixed[[i]]
  XMfix[,-c(index[[i]])] = 0
  Mfix_PC[[i]] = XMfix %*% fixef_PC[[i]] # Matrix of the Groupe effect
  # backtransform the PC to original coefficients
  Mfix[[i]] <- Mfix_PC[[i]]%*%t(spectra_PCA_loadings)
}


######################
# compute the LRT 
######################

# objects initialisation ------------
Res_std_error_PC_null <- vector("list", length=length(null_formulas))
loglik_PC_null <- vector("list", length=length(null_formulas))
sumlog <- c()

### compute the LRT ----------------------------
for (i in 1:length(null_formulas)){
  Res_std_error_PC_null[[i]]  <- sapply(MM_PC_null[[i]], sigma)

  #  sumlog
  loglik_PC_null[[i]] <- sapply(MM_PC_null[[i]], logLik, REML=REML[i])
  
  sumlog[i] <- 2*(sum(loglik_PC_full[i,] - loglik_PC_null[[i]]))
}

names(sumlog) <- names(null_formulas)
sumlog_true <- sumlog

# Graphs ----------------------------

col1="blue"
col2="red"
par(mar=c(4,3,2,6))
dif <- vector(mode = "list")
for (i in 1:length(null_formulas)){
  # graphs
  mat <- cbind(loglik_PC_null[[i]], loglik_PC_full[i,])
  rownames(mat) <- paste0("PC", 1:nPC)
  col <- c(col1,col2)
  par(xpd=TRUE)
  barplot(t(mat), beside=T, ylab="Log-likelihood",
          cex.names=0.8, las=2, col=col, 
          main = paste("log-likelihoods",names(null_formulas)[i]),
          xpd=TRUE)
  legend("topright", legend = c("loglik restricted","loglik full"),
         fill = col, bty = "n", inset=c(-0.1,0))
  
  dif[[i]] <- 2*(mat[,2] - mat[,1])
  
  
}

13.3 Booststrap tests

### Bootstrap ---------------------------------------
# bootstrapLT input arguments:
# MM_null = MM_PC_null$volunteer
# useREML=TRUE
# null_formula <- null_formulas[[i]]

pander("running bootstrap ...")
 
bootstrapLT <- function(useREML, MM_null, null_formula, design, outcomes) {
  
  # simulate y from null models --------
  nPC <- length(MM_null)
  simulatedY <- c()

  for (i in 1:nPC){
    ysim       <- unlist(simulate(MM_null[[i]], re.form=NA))
    simulatedY <- cbind(simulatedY, ysim)
    # dim(simulatedY)
  }

  y <- simulatedY
  dimnames(y) <- dimnames(outcomes)

  # build restricted model -------
  f_null <- parlmer_interaction(design, y, null_formula, REML=useREML)

  # build full model -------
  f_full <- parlmer_interaction(design, y, form, REML=useREML)

  MM_f_null <- f_null$merMod_obj
  MM_f_full <- f_full$merMod_obj

  # LR -------
  loglikelihood_null <- sapply(MM_f_null, logLik, REML=useREML)
 
  loglikelihood_full <- sapply(MM_f_full, logLik, REML=useREML)
 
  ratio <- 2*(loglikelihood_full-loglikelihood_null)
  sumlog <- 2*(sum(loglikelihood_full - loglikelihood_null))
  
  return(list(sumlog=sumlog, ratio=ratio))
  # returns the summed LLR (sumlog) and the LLR per PC (ratio)
}


# test the function bootstrapLT
bootstrapLT(MM_null = MM_PC_null[["Candies" ]], useREML =REML["Candies" ],
                      null_formula = null_formulas[["Candies" ]], 
            design = designInter, outcomes = outcomes)


# Boostrapping: Apply bootstrapLT for each effect
sumlog_boot <- vector("list", length=length(null_formulas))
ratio_boot <- vector("list", length=length(null_formulas))
names(sumlog_boot) <- names(ratio_boot) <- null_effects

set.seed(2018)
for (i in 1:length(null_formulas)){
  null_effect <- null_effects[i]
  res = replicate(nsim, bootstrapLT(MM_null = MM_PC_null[[null_effect]],
               useREML =REML[null_effect],null_formula = null_formulas[[null_effect]],
              design = designInter, outcomes = outcomes),
                          simplify = "array")
  sumlog_boot[[i]] <- res["sumlog", ]
  sumlog_boot[[i]] <- unlist(sumlog_boot[[i]])
  ratio_boot[[i]] <- res["ratio", ]
  ratio_boot[[i]] <- do.call(rbind, ratio_boot[[i]])
}


save(sumlog_boot, sumlog_true,ratio_boot, file=file.path(out_path,name_RData))
load(file=file.path(out_path, name_RData))

13.4 Boostrap p-values

pval <- c()
for (i in 1:length(null_formulas)){
 pval[i] <- (sum(sumlog[i]<sumlog_boot[[i]])+1)/(nsim+1)
}
names(pval) <- names(null_formulas)

pander("p-values")

p-values

pval
##       Candies        Judges CandiesJudges 
##  0.0004997501  0.0354822589  0.0004997501

13.4.1 Graph

difmat <- do.call(cbind, dif)
difmat <- data.frame(PC=substr(rownames(difmat),3,3), difmat)
colnames(difmat) <- c("PC", names(null_formulas))
difmat <- gather(difmat, key=Effect, value = value, Candies , Judges, CandiesJudges)
difmat$Effect <- as.factor(difmat$Effect)
difmat$Effect <- factor(difmat$Effect, levels = c("Candies", "Judges", "CandiesJudges"))

LLRplot <- ggplot(data=difmat, aes(x=PC, y=value, fill = Effect)) +
  geom_bar(width=0.5,stat="identity",
           position=position_dodge(width=0.5)) +
  theme_classic()+labs(title="(Restricted) Log-likelihood Ratios") + 
  ylab(label="(R)LLR") +
  guides(fill=guide_legend(title="Removed effect"))


tab <- paste(c("<", "", "<"), round(pval, 4))
tab <- data.frame(Effect = names(pval), `p-value` = tab, 
                  chi2 = c("<5e-04","-", "-"))
pval_boot.table <- ggtexttable(tab,cols = c("Effect",
                                            "Boostrapped p-value", "Chi2 test"),
                                 rows = NULL)



LLR_pval_plot <- ggarrange(LLRplot, pval_boot.table,
                ncol = 1, nrow = 2,
                heights = c(1,  0.3), common.legend = TRUE)

# ggexport(LLR_pval_plot, filename = file.path(fig_path,"SDA_LLR_pval_plot.pdf"),
#          height = 5, width = 5)

# LLR_pval_plot

# journal article output



p <- LLRplot + annotation_custom(ggplotGrob(pval_boot.table),
                              xmin = 2, ymin = 100)
p

# ggsave(p, filename = file.path(fig_path,"SDA_GLLR.pdf"),
#          width=7, height=4.5, scale=0.9)
### Chapter graph
difmat <- do.call(cbind, dif)
difmat <- data.frame(PC=substr(rownames(difmat),3,3), difmat)
colnames(difmat) <- c("PC", names(null_formulas))
difmat <- gather(difmat, key=Effect, value = value, Candies , Judges, CandiesJudges)
difmat$Effect <- as.factor(difmat$Effect)
difmat$Effect <- factor(difmat$Effect,levels(difmat$Effect)[c(2,1,3)])

difmat2 <- cbind(difmat[difmat$Effect=="Candies", "value"],
                 difmat[difmat$Effect=="Judges", "value"],
                 difmat[difmat$Effect=="CandiesJudges", "value"])
dimnames(difmat2) <- list(c(1:8), levels(difmat$Effect))


# pdf(file.path(fig_path, "SDA_GLLR.pdf"), width=7, 
#     height=6, pointsize = 20)
par(mar=c(4,4,4,4), xpd=TRUE)

# plotting settings -------------------------------------------------------
ylim <- range(mat)*c(1,1.5)
angle1 <- rep(c(45,45,135), length.out=7)
angle2 <- rep(c(45,135,135), length.out=7)
density1 <- seq(5,35,length.out=7)
density2 <- seq(5,35,length.out=7)

op <- par(mar=c(4,3,1,1))


barplot(t(difmat2), beside=TRUE,col = gg_color_hue(3),ylab="(R)LLR", xlab="PC", 
         main = "(Restricted) Log-likelihood Ratios",
        angle=angle1[c(2,4,6)], density=density1[c(2,4,6)], 
        ylim = c(0,200))
barplot(t(difmat2), beside=TRUE, add=TRUE, col = gg_color_hue(3),
        ylab="GLLR", xlab="PC", 
         main = "(Restricted) Log-likelihood Ratios",
        angle=angle2[c(2,4,6)], density=density2[c(2,4,6)], 
        ylim = c(0,200))
legend("topright", legend = c("Candies", "Judges", "Candies*Judges"), 
       title = "Removed effect:", ncol=1, col = gg_color_hue(3),
       fill=gg_color_hue(3),angle=angle1[c(2,4,6)],
       density=density1[c(2,4,6)], inset=c(0,0.1), bty="n")
par(bg="transparent")
legend("topright", legend = c("Candies", "Judges", "Candies*Judges"), 
       title = "Removed effect:",ncol=1, col = gg_color_hue(3),
       fill=gg_color_hue(3),angle=angle2[c(2,4,6)],
       density=density2[c(2,4,6)],inset=c(0,0.1), bty="n")
# dev.off()


# pdf(file.path(fig_path, "SDA_pval.pdf"))
pval_boot.table
# dev.off()

13.5 Histogram drawings

names(sumlog_boot) <- casefold(names(null_formulas), upper = FALSE)
df <- nPC *   4    
######################
# Candies
######################
# pdf(file = file.path(fig.path, "hist_time.pdf"), width = 7, height = 5)
par(mar=c(2.1, 2.1, 2.1, 7.5), xpd=TRUE, mfrow=c(1,1))

m=hist(sumlog_boot$candies, freq=F, breaks=100, 
                 xlab="Global Likelihood Ratio Statistic",
                 xlim=range(sumlog["Candies"], sumlog_boot$candies),
                 ylim = c(0,0.08),
                 col = gray67,border = gray67, 
                 main = " Fixed Candies effect", cex.main = 2.2)
            lines(density(sumlog_boot$candies), col=darkblue,lwd=2)
            lines(dchisq(seq(0,max(sumlog_boot$candies)), df),
                  col= "limegreen", lwd = 4)
            points(sumlog["Candies"], 0, col="red", pch=19, lwd=6)
            legend("topright", 
                   legend = c(paste0("True GLLR: ", round(sumlog["Candies"],2)), 
                              "Kernel density", 
                              paste0("chi2 distrib. (df=", df,")")), 
                   col = c("red", darkblue, "limegreen"),
                   lty=c(NA,1,1),pch=c(19,NA,NA),
                   inset=c(-0.2,0),box.lty=0, cex = 1.4, 
                   y.intersp = 0.8, lwd=c(4))
            
# dev.off()      
# 

######################
# Judges
######################

# pdf(file = file.path(fig.path, "hist_volunteer.pdf"), width = 7, height = 5)
# df*(1/2*dchisq(0:30, 1)+ (1/2*dchisq(0:30, 0)))
            
par(mar=c(2.1, 2.1, 2.1, 7.5), xpd=TRUE, mfrow=c(2,1))
df <- nPC 
m=hist(sumlog_boot$judges, freq=F, breaks=100, 
                 xlab="Global Likelihood Ratio Statistic",
                 xlim=range(sumlog["Judges"], sumlog_boot$judges),
                 col = gray67,border = gray67, 
                 main = " Random Judges effect", cex.main = 2.2)
            lines(density(sumlog_boot$judges), col=darkblue,lwd=2)
            # lines(0:30,(1/2*dchisq(seq(0,30), df)+ (1/2*dchisq(seq(0,30), 0))),
            #        col= "limegreen", lwd = 4)
            points(sumlog["Judges"], 0, col="red", pch=19, lwd=6)
            legend("topright", 
                   legend = c(paste0("True GRLLR: ", round(sumlog["Judges"],2)), 
                              "Kernel density"), 
                   col = c("red", darkblue), lty=c(NA,1),pch=c(19,NA),
                   inset=c(-0.2,0),box.lty=0, cex = 1.4, 
                   y.intersp = 0.8, lwd=c(4))
# dev.off()  

######################
# CandiesJudges
######################

# pdf(file = file.path(fig.path, "hist_sampling.pdf"), width = 7, height = 5)
par(mar=c(2.1, 2.1, 2.1, 7.5), xpd=TRUE, mfrow=c(1,1))

m=hist(sumlog_boot$candiesjudges, freq=F, breaks=100, 
                 xlab="Global Likelihood Ratio Statistic",
                 xlim=range(sumlog["CandiesJudges"], sumlog_boot$candiesjudges),
                 col = gray67,border = gray67, 
                 main = " Random CandiesJudges effect", cex.main = 2.2)
            lines(density(sumlog_boot$candiesjudges), col=darkblue,lwd=2)
            points(sumlog["CandiesJudges"], 0, col="red", pch=19, lwd=6)
            legend("topright", 
                   legend = c(paste0("True GRLLR: ",
                                     round(sumlog["CandiesJudges"],2)), 
                              "Kernel density"), 
                   col = c("red", darkblue), lty=c(NA,1),pch=c(19,NA),
                   inset=c(-0.2,0),box.lty=0, cex = 1.4, 
                   y.intersp = 0.8, lwd=c(4))
# dev.off()      
# 

names(sumlog_boot) <- casefold(names(null_formulas), upper = FALSE)
df <- nPC *   4    
######################
# Candies
######################
# pdf(file = file.path(fig_path, "SDA_hist_Candies.pdf"),
#     width = 7, height = 5, pointsize = 15)

par(mar=c(2.1, 2.1, 2.1, 3), xpd=TRUE, mfrow=c(1,1))

m=hist(sumlog_boot$candies, freq=F, breaks=100, 
                 xlab="Global Likelihood Ratio Statistic",
                 xlim=range(pretty(c(sumlog["Candies"],
                            sumlog_boot$candies))),
                 ylim = c(0,0.06),
                 col = "gray75",border = "gray75", 
                 main = "Candies", cex.main = 2.2)
            lines(density(sumlog_boot$candies), col=darkblue,lwd=2,
                  lty=2)
            lines(dchisq(seq(0,max(sumlog_boot$candies)), df), 
                  col= "chartreuse4", lwd = 2)
            points(sumlog["Candies"], 0, col="red", pch=19, lwd=6)
            legend("topright", 
            legend = c(paste0("True GLRT: ", 
                          round(sumlog["Candies"],2)), 
                              "Kernel density", 
                              paste0("chi2 distrib. (df=", df,")")), 
                   col = c("red", darkblue, "chartreuse4"),
                   lty=c(NA,2,1),pch=c(19,NA,NA),
                   inset=c(-0.1,0),box.lty=0, cex = 1.4, 
                   y.intersp = 0.8, lwd=c(4))
            
# dev.off()

######################
# Judges
######################

# pdf(file = file.path(fig_path, "SDA_hist_Judges.pdf"),
# width = 7, height = 5, pointsize = 15)

# df*(1/2*dchisq(0:30, 1)+ (1/2*dchisq(0:30, 0)))
            
par(mar=c(2.1, 2.1, 2.1, 4.5), xpd=TRUE)
df <- nPC 
m=hist(sumlog_boot$judges, freq=F, breaks=100, 
                 xlab="Global Likelihood Ratio Statistic",
                 xlim=range(sumlog["Judges"], sumlog_boot$judges),
                 col = "gray75",border = "gray75", 
                 main = "Judges", 
       cex.main = 2.2)
            lines(density(sumlog_boot$judges), col=darkblue,lwd=2,
                  lty=2)
            # y <- (1/2*dchisq(seq(0,30), df)+ 
            #         (1/2*dchisq(seq(0,30), 0)))
           
            # lines(0:30,y, col= "limegreen", lwd = 4)
            
            points(sumlog["Judges"], 0, col="red", pch=19, lwd=6)
            legend("topright", 
                   legend = c(paste0("True GLRT: ", round(sumlog["Judges"],2)), 
                              "Kernel density"), 
                   col = c("red", darkblue), lty=c(NA,2),pch=c(19,NA),
                   inset=c(-0.2,0),box.lty=0, cex = 1.4, 
                   y.intersp = 0.8, lwd=c(4))
# dev.off()

######################
# CandiesJudges
######################

# pdf(file = file.path(fig_path, "SDA_hist_CandiesJudges.pdf"), width = 7, height = 5, pointsize = 15)
par(mar=c(2.1, 2.1, 2.1, 3), xpd=TRUE, mfrow=c(1,1))

m=hist(sumlog_boot$candiesjudges, freq=F, breaks=100, 
                 xlab="Global Likelihood Ratio Statistic",
                 xlim=range(pretty(c(sumlog["CandiesJudges"],
                            sumlog_boot$candiesjudges))),
                 col = "gray75",border = "gray75", 
                 main = "Candies*Judges", cex.main = 2.2)
            lines(density(sumlog_boot$candiesjudges), col=darkblue,lwd=2,
                  lty=2)
            points(sumlog["CandiesJudges"], 0, col="red", pch=19, lwd=6)
            legend("topright", 
                   legend = c(paste0("True GLRT: ",
                                     round(sumlog["CandiesJudges"],2)), 
                              "Kernel density"), 
                   col = c("red", darkblue), lty=c(NA,2),pch=c(19,NA),
                   inset=c(-0.1,0),box.lty=0, cex = 1.4, 
                   y.intersp = 0.8, lwd=c(4))
# dev.off()

13.6 Compare the theorerical and true quantiles for fixed effect

13.6.1 Graph histo per PC

# pdf(file = file.path(fig_path, "SDA_hist_perPC_Candies.pdf"),
# width = 11, height = 10, pointsize = 10)

par(mar=c(2, 4, 4, 2), xpd=TRUE, mfrow=c(3,3), cex=1)
ratio_boot_Candies <- ratio_boot$Candies
for (i in 1:dim(ratio_boot_Candies)[2]){
  m=hist(ratio_boot_Candies[,i],  freq=F, breaks=100, 
                 xlab=paste0("LLR"),
                 col = "gray75",border = "gray75",
                 main = paste0("LLR Candies - PC ",i))
            lines(density(ratio_boot_Candies[,i]), col=darkblue,lwd=3,
                  lty=2)
            curve(dchisq(x, df=4), col='chartreuse4', 
                   main = "Chi-Square Density Graph", 
                  from=0,to=30, add=TRUE, lwd=3, xpd=F)
            
            # points(sumlog["CandiesJudges"], 0, col="red", pch=19, lwd=6)
            
     
}
plot(NULL, xlim=c(0,1), ylim=c(0,1), bty="n", xaxt="n", yaxt="n",
     xlab="", ylab="")
legend("topright", legend = c("Kernel density", "chi2 (df=4)"),
                    col = c(darkblue,  "chartreuse4"), lty=c(2,1),
                    inset=c(0,0),box.lty=0,
                   y.intersp = 1, lwd=c(4), cex=1.4)

# dev.off()

13.6.2 Compare the quantiles by PC and for the GLLR

# By PC
df <- 4 
theo_quant <- qchisq(seq(0, 0.9999, 1/nsim), df=df)
# pdf(file = file.path(fig_path,
#      paste0("SDA_RLLRqqplot_perPC_","Candies",".pdf")),
#       width = 12, height = 6)
  par(mfrow=c(2,4), cex=1, mar=c(4,4,4,1))
  for (i in 1:dim(ratio_boot_Candies)[2]){
    boot_quant <- quantile(ratio_boot_Candies[,i],probs = seq(0, 0.9999, 1/nsim))
    qTest <- quantile(ratio_boot_Candies[,i],probs = seq(0.0001,0.999, 1/nsim))
    
    qChi2 <- theo_quant
    # boot_quant <- subset(boot_quant, !names(boot_quant) %in% c("0.00%", "99.95%","100.00%"))
    plot(boot_quant, qChi2,cex.main=1,
         xlim=range(boot_quant),
         ylim=range(qChi2),
         main = paste("LLR Candies -", colnames(ratio_boot_Candies)[i]), xlab="Sample quantiles",
         ylab="theoretical quantiles", cex.main=1)
    lines(c(0, min(max(boot_quant), max(qChi2))),
          c(0, min(max(boot_quant), max(qChi2))),
          col="red", lwd=1.5)
  }

# dev.off()

# Globally
df <- nPC*4 
Theo_quant <- qchisq(seq(0, 0.9999, 1/nsim), df=df)
boot_quant <- quantile(sumlog_boot$candies,probs = seq(0, 0.9999, 1/nsim))
par(mfrow=c(1,1))

# pdf(file = file.path(fig_path, "SDA_GLLRquantiles_Candies.pdf"), 
#     width = 6, height = 6, pointsize = 15)

plot(boot_quant,Theo_quant, main = "True and theoretical GLLR for Candies effect", xlab = "Sample quantiles", ylab = "Theoretical quantiles", xlim = range(boot_quant, Theo_quant),
     ylim=range(boot_quant, Theo_quant))
lines(x=c(1:max(boot_quant, Theo_quant)), y=c(1:max(boot_quant, Theo_quant)), col="red")

# dev.off()

13.6.3 p-value based on chi2 distribution

pander("true GLLR:")

true GLLR:

sumlog["Candies"]
##  Candies 
## 284.9536
curve(dchisq(x, df=4*nPC), col='red', main = "Chi-Square Density Graph",
          from=0,to=60)

pander("p-value")

p-value

pchisq(sumlog["Candies"], df=4*nPC, lower.tail=FALSE)
##      Candies 
## 2.294577e-42

13.7 Compare the theorerical and true quantiles for random effects

13.7.1 Graph histo per PC

# pdf(file = file.path(fig_path, "SDA_hist_perPC_Judges.pdf"), 
#     width = 10, height = 10, pointsize = 10)
par(mar=c(2, 4, 4, 3),  mfrow=c(3,3), cex=0.9)

ratio_boot_Judges <- ratio_boot$Judges
for (i in 1:dim(ratio_boot_Judges)[2]){

x <- seq(0,round(max(ratio_boot_Judges[,i])),0.001)

y <- 0.5*(dchisq(x, df=1) + dchisq(x, df=0))
yy <- dchisq(x, df=1)

  m=hist(ratio_boot_Judges[,i],  freq=F, breaks=100, 
                 xlab=paste0("RLLR"),
                 col = "gray75",border = "gray75",
                 main = paste0("RLLR Judges - PC ",i))
 
 lines(density(ratio_boot_Judges[,i]), col=darkblue,lwd=2,
       ylim =c(0,7), lty=2)
 
 lines(x,y, type="l", col="red", lwd=2,ylim =c(0,7),xpd=F, lty=3)
 
 lines(x,yy, type="l", col="chartreuse4", lwd=2,ylim = c(0,7),xpd=F)
 

     
}
plot(NULL, xlim=c(0,1), ylim=c(0,1), bty="n", xaxt="n", yaxt="n",
     xlab="", ylab="")
  legend("topright",legend = c("Kernel density", 
                               "mixture of chi2\n(df=0,1)", 
                               "chi2 (df=1)"), 
     col = c(darkblue,  "chartreuse4", "red"), lty=c(2,1,3),
     inset=c(0,0),box.lty=0, 
     y.intersp = 0.8, lwd=c(4),xpd=TRUE)

# dev.off()

par(mfrow=c(1,1))

# pdf(file = file.path(fig_path, "SDA_hist_perPC_CJ.pdf"), 
#     width = 10, height = 10, pointsize = 10)

par(mar=c(2, 4, 4, 3), mfrow=c(3,3), cex=0.9)

ratio_boot_CJ <- ratio_boot$CandiesJudges
for (i in 1:dim(ratio_boot_CJ)[2]){

x <- seq(0,round(max(ratio_boot_CJ[,i])),0.001)

y <- 0.5*(dchisq(x, df=1) + dchisq(x, df=0))
yy <- dchisq(x, df=1)

  m=hist(ratio_boot_CJ[,i],  freq=F, breaks=100, 
                 xlab=paste0("RLLR"),
                 col = "gray75",border = "gray75",
                 main = paste0("RLLR for CJ - PC ",i))
 
 lines(density(ratio_boot_CJ[,i]), col=darkblue,lwd=2, lty=2)
 
 lines(x,y, type="l", col="red", lwd=2,xpd=F, lty=3)
 
 lines(x,yy, type="l", col="chartreuse4", lwd=2,xpd=F)
  

     
}

plot(NULL, xlim=c(0,1), ylim=c(0,1), bty="n", xaxt="n", yaxt="n",
     xlab="", ylab="")
  legend("topright",legend = c("Kernel density", 
                               "mixture of\nchi2 (df=0,1)", 
                               "chi2 (df=1)"), 
     col = c(darkblue,  "chartreuse4", "red"), lty=c(2,1,3),
     inset=c(0,0),box.lty=0, 
     y.intersp = 0.8, lwd=c(4),xpd=TRUE)

# dev.off()  

13.7.2 Compare the quantiles

13.7.2.1 Compute theo quant

f <- function(x, P){
  0.5*pchisq(x,0)+0.5*pchisq(x,1)-P
}
qchisq(0.1,0)
## [1] 0
vect <- c()
for (i in seq(0.0001,0.999,by=1/nsim)){
  # print(i)
  interval <- c(min(qchisq(i,0),qchisq(i,1)), 
                max(qchisq(i,0),qchisq(i,1)))
  interval <- c(0,10)
  vect <- c(vect, uniroot(f,interval, 
                          tol = 0.0001, P = i)$root)
}
theo_quant <- vect
# plot(theo_quant)

### Per PC
############
# bootval <- ratio_boot[["Judges"]]
# bootval <- ratio_boot[["Assessors"]]
# bootval <- ratio_boot[["CandiesJudges"]]

names(ratio_boot) <- c( "Candies","Judges" , "CJ")
 
for (j in 2:3){# random effects
  # pdf(file = file.path(fig_path,
  #    paste0("SDA_RLLRqqplot_perPC_",names(ratio_boot)[j],".pdf")),
  #     width = 12, height = 6)

    par(mfrow=c(2,4), cex=1, mar=c(4,4,4,1))
  for (i in 1:dim(ratio_boot[[j]])[2]){
    qTest <- quantile(ratio_boot[[j]][,i],probs = seq(0.0001,0.999, 1/nsim))
    qChi2 <- theo_quant
    qTest <- subset(qTest, !names(qTest) %in% c("0.00%", "99.95%","100.00%"))
    plot(qTest, qChi2, 
         xlim=range(qTest),
         ylim=range(qChi2),
         main = paste("RLLR for", names(ratio_boot)[j], "-", 
                      colnames(ratio_boot[[j]])[i]), xlab="Sample quantiles",
         ylab="theoretical quantiles", cex.main=1)
    lines(c(0, min(max(qTest), max(qChi2))),
          c(0, min(max(qTest), max(qChi2))),
          col="red", lwd=1.5)
  }
  
  # dev.off()
}

14 Session info

sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] doBy_4.6.6          car_3.0-8           carData_3.0-4      
##  [4] mdatools_0.10.3     tidyr_1.1.0         ggpubr_0.3.0       
##  [7] emdbook_1.3.12      spatstat_1.64-1     rpart_4.1-15       
## [10] nlme_3.1-148        spatstat.data_1.4-3 reshape2_1.4.4     
## [13] cowplot_1.0.0       ggplot2_3.3.2       gridExtra_2.3      
## [16] stringr_1.4.0       pander_0.6.3        plyr_1.8.6         
## [19] lme4_1.1-23         Matrix_1.2-18       MBXUCL_0.1         
## [22] here_0.1           
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1            backports_1.1.8         splines_4.0.1          
##   [4] crosstalk_1.1.0.1       digest_0.6.25           foreach_1.5.0          
##   [7] htmltools_0.5.0         magrittr_1.5            tensor_1.5             
##  [10] cluster_2.1.0           openxlsx_4.1.5          recipes_0.1.13         
##  [13] gower_0.2.2             bdsmatrix_1.3-4         colorspace_1.4-1       
##  [16] haven_2.3.1             xfun_0.15               dplyr_1.0.0            
##  [19] crayon_1.3.4            jsonlite_1.7.0          survival_3.2-3         
##  [22] iterators_1.0.12        ape_5.4                 glue_1.4.1             
##  [25] polyclip_1.10-0         gtable_0.3.0            ipred_0.9-9            
##  [28] R2HTML_2.3.2            webshot_0.5.2           spls_2.2-3             
##  [31] BiocGenerics_0.34.0     abind_1.4-5             scales_1.1.1           
##  [34] mvtnorm_1.1-1           rstatix_0.5.0           miniUI_0.1.1.1         
##  [37] Rcpp_1.0.4.6            xtable_1.8-4            ropls_1.20.0           
##  [40] foreign_0.8-80          proxy_0.4-24            stats4_4.0.1           
##  [43] lava_1.6.7              prodlim_2019.11.13      htmlwidgets_1.5.1      
##  [46] ellipsis_0.3.1          pkgconfig_2.0.3         farver_2.0.3           
##  [49] nnet_7.3-14             deldir_0.1-25           caret_6.0-86           
##  [52] tidyselect_1.1.0        labeling_0.3            rlang_0.4.6            
##  [55] manipulateWidget_0.10.1 later_1.1.0.1           munsell_0.5.0          
##  [58] phyclust_0.1-29         cellranger_1.1.0        tools_4.0.1            
##  [61] generics_0.0.2          ade4_1.7-15             pls_2.7-2              
##  [64] broom_0.5.6             evaluate_0.14           fastmap_1.0.1          
##  [67] yaml_2.2.1              goftest_1.2-2           ModelMetrics_1.2.2.2   
##  [70] knitr_1.29              zip_2.0.4               rgl_0.100.54           
##  [73] purrr_0.3.4             mime_0.9                compiler_4.0.1         
##  [76] curl_4.3                e1071_1.7-3             ggsignif_0.6.0         
##  [79] spatstat.utils_1.17-0   tibble_3.0.1            statmod_1.4.34         
##  [82] stringi_1.4.6           forcats_0.5.0           lattice_0.20-41        
##  [85] nloptr_1.2.2.1          vctrs_0.3.1             pillar_1.4.4           
##  [88] lifecycle_0.2.0         data.table_1.12.8       httpuv_1.5.4           
##  [91] R6_2.4.1                promises_1.1.1          clusterSim_0.48-3      
##  [94] rio_0.5.16              codetools_0.2-16        boot_1.3-25            
##  [97] MASS_7.3-51.6           rprojroot_1.3-2         withr_2.2.0            
## [100] Deriv_4.0               mgcv_1.8-31             parallel_4.0.1         
## [103] hms_0.5.3               timeDate_3043.102       coda_0.19-3            
## [106] class_7.3-17            minqa_1.2.4             clValid_0.6-8          
## [109] rmarkdown_2.3           bbmle_1.0.23.1          pROC_1.16.2            
## [112] numDeriv_2016.8-1.1     Biobase_2.48.0          shiny_1.5.0            
## [115] lubridate_1.7.9